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Abstract 

We present extensive Monte Carlo simulations for the thermodynamic and structural properties 
of a planar bilayer of dipolar hard spheres for a wide range of densities, dipole moments and layer 
separations. Expressions for the stress and pressure tensors of the bilayer system are derived. 
For all thermodynamic states considered the interlayer energy is shown to be attractive and 
much smaller than the intralayer contribution to the energy. It vanishes at layer separations of 
the order of two hard sphere diameters. The normal pressure is negative and decays as a function 
of layer separation /i as — Intralayer and interlayer pair distribution functions and angular 
*author for correspondance : martial.mazars@th.u-psud.fr 
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correlation functions arc presented. Despite the weak interlayer energy strong positional and 
orientational correlations exist between particles in the two layers. 
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I INTRODUCTION 

Dipolar interactions play a significant role in determining the structural, magnetic or rheo- 
logical properties of a variety of quasi two-dimensional (2D) systems (monolayers, multilayers, 
thin films) including suspensions of colloidal particles at an air-water interface, adsorbed am- 
phiphilic molecules, lipid bilayers, ultrathin magnetic films etc.. (see e. g. ref. [1 and references 
therein). In most of these systems the properties and phase behavior result, though, from an 
interplay of the dipolar interaction with competing interactions, as for instance, the hydrocar- 
bon chain tails or water mediated interactions in lipid bilayers [21 E], or exchange interaction 
and magneto crystalline anisotropy in thin magnetic films [4J. Although simulations taking into 
account full atomic details have been performed in the past (generally computationally costly) 
for these kinds of systems (see e.g. ref. [5] and references therein) we believe that a study of a 
purely dipolar bilayer system is of interest in its own right providing unbiased insight into the 
role of the dipolar interaction. The experimental system which perhaps comes closest to the pure 
dipolar system is the ferrofluid system. In effect, association into chains, rings, branched struc- 
tures or stripes has been demonstrated in recent experiments on strongly interacting (Fe304) 
ferrofluids [H 13 [8] and comparison with simulation results presenting similar structures is more 
than suggestive that the dipolar hard sphere (DHS) system is a fair representation of these types 
of ferrofluid. 

Extensive Monte Carlo (MC) simulation and theoretical results for the self organization of quasi 
2D DHS are already available for the monolayer system both with and without an external field 

iniiiniiiiiiiziiiaiiaiiaiis]. 

The purpose of the present paper is to extent these results to a symmetric planar bilayer the 
main interest, evidently, being to probe the effect of the interlayer interaction on particle orga- 
nization. 

In Sect, [n] we define the bilayer model and give details of the numerical simulation methods we 
use. The next section gives expressions for the energy, stress tensor and correlation functions of 



the bilayer system. Sect. IV contains the simulation results for the thermodynamic and struc- 
tural properties. A summary is given in the last section. The three appendices A-C provide 
expressions for the Ewald sums of energy (A), pressure and forces (C) and a derivation of the 
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microscopic stress tensor of the bilayer (B). 



II MODEL AND NUMERICAL METHODS 

The systems consist of = 2A''o particles with permanent point dipole moment fi interacting 
via hard sphere and dipolar potentials. Particles are evenly distributed among two layers Li 
and L2 separated by a distance h, each layer being rectangular with sides and Ly ; A = L^Ly 
is the surface area of the layers. Periodic boundary conditions (p.b.c), with spatial periodicities 
Lx and Ly, are applied in the directions x and y parallel to the layers, but no p.b.c. are taken in 
the third direction z. Particle positions are constrained to lie in the layers but dipole moments 
can orient in full 3D space. The interaction potential between the particles is pairwise additive 
and is represented as 



00 for Tij < a 

(1) 
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Hi ■ - 3{Hi ■ rij){fij ■ rij)\ for nj > a 



where a = 1 is the hard sphere diameter taken as unit length, /Xj the dipole moment of particle 
i and fjj = rij/rij the unit bond vector between particles i and j. In the following, we will use 
the notations 

Tij = Sij + Zij and fi^ = ^/x^ (2) 

where is the unit vector perpendicular to the layers and /xj a unit vector in the direction of 
dipole moment i. 

Only surface separation h > 1 which avoid hard core interactions between the layers have been 
considered. A few simulation results for h > 1 have been presented previously by one of us [1]. 
Monte Carlo (MC) simulations have been performed in the canonical (NVT) ensemble with 
system sizes comprising N = 1024 — 3200 particles. The total number of MC cycles varied 
from 0.2 x 10^ to 2 x 10^, depending on density and dipole moment, each cycle consisting of 
displacement and rotation of the N particles. The amplitude of the trial moves was chosen to 
obtain acceptance ratios between 30 and 50% for each thermodynamic state. No exchange of 
particles between layers Li and L2 is allowed. 

Reduced quantities for surface area. A* = Ajo^ , surface density p* = pcr^ = N^jAa'^, and dipole 



moment fi* = {^'^ /kTa^Y^"^ will be used throughout the paper. For notational convenience the 
stars will be dropped. 



Ill THERMODYNAMICAL AND STRUCTURAL 

QUANTITIES 

A Energy 

In our model the energy of the bilayer is entirely given by the dipolar contribution which we 
split into an intralayer contribution, JJ^"-^^'"-^ and an interlayer contribution, f/™*'^''^ as 

Udd = U'''^'''' + C/*"*"^ (3) 

These are computed using the Ewald method |17| [TBI [1] [19] ; the relevant expressions for JJ^^^^"- 
and f/*"*^'' are given in Appendix A. 

For bulk systems with slab geometry where periodicity applies only in two spatial directions, say 
Lx and Ly, the Ewald sums are computationally costly due to the appearance in the reciprocal 
space term of a double sum over the distance Zij in the bounded direction of particles i and 
j \17\ [T8] . As in the present case the distance Zij between two particles will be constant, the 
corresponding sums can be reduced to order similar to the cases of Coulomb |20| |2T] or 

Yukawa [22] potentials. 

One can note that the 3D bilayer system can be mapped onto a two-component monolayer 
system by considering the particles in the two layers as distinct species |23]. For most of the 
thermodynamical and structural quantities, both approaches are equivalent ; for instance, in 
the two-component monolayer, is the total interaction between particles belonging to 

different species (different layers). As outlined in the next subsection and in Appendix B, for 
pressures and stresses such a mapping is slightly less straightforward. 

B Surface stress tensor and normal pressure 

Characterizing the pressure in the bilayer system needs some care. In particular, since 
the particles are constrained to belong to layers Li and L2, some degrees of freedom of the 
particles are frozen by the geometrical features of the system. These constraints have obviously 



an influence on the flux of momentum per unit area in the system and therefore affect the stress 
tensor. For the sake of definitness a full derivation of the stress tensor from the lagrangian 
function of the bilayer system is given in Appendix B. 

As for systems with slab geometry or interfaces ||24i , the stress tensor is decomposed into lateral 



and normal components. According to Eq.(B.13 B.15), the lateral component to the pressure 
tensor is given by 



(4) 



E E ^ir^M^ij^h)'^ 
where ^{sij, h) is the pair potential. 

From the point of view of mapping the bilayer system onto a two-component monolayer system, 
the lateral pressure Ut in the bilayer, defined in Eq.Q through Eqs.(B.12 B.15 ), corresponds to 
the pressure of the 2D, two-component monolayer system. In solid surface physics, Ut is related 



to the surface stress fj by 11^ = —77 (cf. Eq.(B.15)), and for fluids confined in slab geometry Ut 
is related to the lateral pressure Pt{z) by 

Ut= dz Pt{z). 



Ut can be composed into ideal, hard sphere (HS), and dipolar contributions 

Ut = 2pkT + 2n^^^) + ng2^^ + n^''^) 

where the dipolar part 11^'^^ is obtained from Eq.(jlJ and the relation 



(5) 



^V2^(<^\s,„h) =3 



{sl + /l2)5/2 



Mi /^, -5- 



4 + ^' 



(6) 
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(4 + /l2)5/2 



(see Eq. (B.14) of Appendix B). n^'^^ contains both intralayer contributions of layers Li and 



L2 and the interlayer contribution; thus, for h — > cxd, 11^°'^ is twice the dipolar contribution to 
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the 2D pressure of a monolayer. The dipolar inter layer contribution to Ht is given by the last 
contribution in the right-hand side (r.h.s.) of Eq.Q ; this contribution becomes very small as 
soon as h>2. 

The hard sphere contributions n^"^^ and 11^^]^^, are computed from the contact values of the 
intralayer, gintrai'^)^ ^^"^ interlayer, ginter^ P^i^ distribution functions, defined below, as 



4^^) = "^p'kTg^^S^M 



(7) 



As in the present work, /i > 1 in all computations, we always have H^^^^^^ = 0. In the limit 
h ^ oo and ^ — > 0, 11^ equals the excess contribution to the pressure of a monolayer of hard 
disks with surface density p. Moreover, for h > 1 and /x = 0, 11^ can be approximated quite 
accurately by available equations of state of hard disks (see e.g. ref . [25] ) . 

The asymptotic behaviour of IIt given by Eq.(|5]) can be understood as follows. In the limit 
/i — > oo and fi 0, Ht, given by Eq.([5]), is exactly twice the 2D pressure of a monolayer 
of DHS with the same p and p. In this limit, if the system is viewed as a two component 
monolayer system, the two species remain distinct but there will be no interaction between 
particles belonging to different species. Thus, 117-/2 is exactly the partial pressure of each 
component and the bilayer is fully equivalent to a mixture of two kinds of particles confined in a 
monolayer with HS and dipolar interactions between like particles but no interactions between 
unlike particles. 

In the opposite limit h ^ and p ^ 0, the two species become equal and the bilayer system 
reduces to a one component monolayer system with a surface density 2p (provided that 2p is 
less than the density at close packing of hard disks). Obviously, in this limit, the contribution 
inter included in Eq.(5), and 11-^ equals the 2D pressure of a monolayer of 

dipolar hard disks with a surface density 2p and same p. Also, as in this limit particles become 
indistinguishable, entropy contributions must be modified accordingly. 
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The average normal force by unit area (or normal pressure) is obtained from Eq.(B.19) as 

=-i(a^EE-^(%>^)U) = -:4( dh ) 



p{dd) , p(ffS) 



where Pjf'^^ and Pz^^^ denote the contributions from dipolar and HS interactions, respectively. 
The dipolar parts, Pzt'^^ and n^'^\ are computed using Ewald sums, as described in Appendix 
C. Since in the present work all computations are done with h > 1 one has always Pzz = 0. 
The HS repulsion does, however, contribute to the normal component of the pressure tensor 
indirectly via the spatial positions of the particles in the layers. A similar remark applies to 
the interlayer correlation functions defined below. Eq.([8]) agrees with previous derivations for 
the normal pressure in slab-like geometry |26| [271 EE] or interfaces [23] . The main difference 
between Eq.Q and these relations is that there is no kinetic (ideal gas) contribution in Eq.([8]), 



as a consequence of the constraints that apply to the bilayer systems (cf. Eq.(B.7)). Thus, Pzz 
has to be considered as an average force by unit area normal to the surface rather than a normal 
pressure. 

The surface stress tensor is related to the surface free energy par unit area 7 (or surface tension) 
by the Shuttleworth equation [29] 

rjalB = l^alS + (9) 



where is the 2D strain tensor. In fluid phases, the second contribution in the r.h.s. of Eq.(|9| 
is null and Eq.(|9]) reduces to rjafs = 7^0/3. This is the case in most computations done in the 
present work, except those at high densities. Since in our computations the surface and the 
shape of the layers are kept constant, we do not have access to 7. 



C Correlation functions 

The structure of the bilayer system has been characterized, analogously to the monolayer 
case [141 IS] > by a one particle orientational distribution function of the dipoles and several pair 
correlation functions. 
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The orientational distribution function f{fi), measuring the orientation of the particle dipole 
moments with respect to the layer normal, is defined from the one-body density as 

P 



pW(r, A) = - r)<5(A. - A)) = 



(10) 



(11) 



Pair correlation functions are derived from the general definition of the two-body density 

N 

(r, r', A, AO = ( E <^(^^ - - - mi^j - A') 

where A A' unit vectors along the dipole moments. Specifying to intralayer and 

(2) 

interlayer Plj^^^ two-body surface densities one has 



i&Li j&L-i_,j^i 



+ ^ XI I D'^^Aj - A)<^(Aj - A'; 

ieL2 3&L2,j^i 



(12) 



pSSer(«' A, A') = ^( X X I D'^^Ai - A)^(Ai - A')^ 

(12) distribution functions are related to the two- 



The intralayer 5mtra( 12) and interlayer ginter 
body densities through 

5mtm(12) = 1 + hintra{^2) = \ — ) Pintrai^^ Al, A2) 
5(inter(12) = 1 + hinter{l2) = [ — ) Pinteri^^ Al, A2) 

In particular, the intralayer gf^trai^) ^^'^ interlayer ginteri^) center-to-center pair distribution 
functions are given by 



(13) 



iGLi j&LiJj^i i&L2j&L2,j^i 



Sij 



dintra (12) 



A1A2 



(14) 



ieLi jeL2 



where Sj is the in-plane position of particle i according to the notations defined in Eq.([2]) and 
< ■ >#"i n denotes averaging over orientations of the dipole moments. The angular dependent 
pair correlation functions /i(12) have been expanded, as usual, on a basis set of rotational 
invariants i'^'^i |3ol [31] 

^(12) = Hh, I2, 1; r)l>'i'^'(Ai, A2, ^) (15) 

ll,l2,l 

where the ^'i'^' g^j-g related to the standard rotational invariants ^'i'^/ jj^ an expansion on 
spherical harmonics by (see e.g. [32] ) 

^- \o y 

The most significant projections of the intralayer hintrai^'^) and interlayer hinteri^'^) correlation 
functions calculated in this work are those onto ^^^^ , $112 g^j^^ ^220_ ^]^q correponding expres- 
sions are summarized in Table HI 



D Order parameter 

Possible orientational (nematic) order in a layer can be established from the non-vanishing 
of the second-rank order parameter P2 calculated as the average value of the largest eigenvalue 
of the matrix [HS] 

1 ^° 1 

Q"f^ = jr^2 ^^^"^^ " '^"^^ ' ^^^^ 

where /x^ is the a component of the unit vector /Xj- One can note that the projection h'^^^ obeys 
the asymptotic relationship 

/i220(s) ^ 5p2^ s ^ 00 (18) 
As will be shown below no global nematic order occurs in the systems for p < 0.7. 

IV RESULTS 

A One-body orientational distribution function 

One-body distribution functions /(/i) = /(cos(^)), with polar angle 6 defined by cos6 = 
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fi ■ Bz, obtained from MC simulation at various thermodynamic states are shown in Fig. [TJ^a). 
It is seen that for all states an excellent fit to the MC data is obtained with the one parameter 
function 

/(cos 0; a) = /o exp(— a cos^ 9) (19) 

with normalization constant 

/o = xl^^r-^ (20) 
Values of a obtained by fitting the MC histograms P{cos9), normalized to one, are given in 
Tables [IT IV The results for the orientational distribution functions of the bilayer system are 



quite similar to those obtained earlier for monolayers [T2]. As /x increases the dipole moments 
tilt more and more into the layer plane {cos 9 ~ 0). The interaction between the two layers 
induces, though, a slight effect, in comparison to the monolayer system, as seen in Fig. [l]^b) 
showing the variation of the orientational distributions with inter layer separation h for p = 0.7 
and /U = 2.00. As the separation between the layers decreases, the coupling between layers 
increases which entails a slight tendency of the dipoles to orient perpendicularly to the plane. 
As a consequence the distributions are slightly broadened (the value of a decreases). 

B Energy 

The variation of the intralayer PW^^'^"' / N and interlayer /3C/*"'**^^/A'" energies as a function of 
layer separation are summarized in Table [IT] for the density p = 0.7 and the two dipole moments 
p = 1 and 2. The intralayer energy is seen to be by far the dominant contribution and is 
nearly independent of h especially at the largest dipole moments where in-plane orientation of 
the dipole moments is prevalent. The interlayer energy is much smaller and decreases rapidly 
with layer separation vanishing at /i ~ 2. The total energy remains practically constant when h 
varies from 1.05 to 2.0. 

Attard and Mitchell have applied a second order perturbation theory on a bilayer of orientable 
dipoles \33\ [M] and found that the interaction free energy between the surfaces decays as the 
fourth power of h at large separation. An analysis of our MC data, for h > 1.6, agrees with 
the behavior obtained in the computations done by Attard and Mitchell ; more precisely, the 
variation of the interlayer energy with h, for p = 0.7 and p=l and 2, can be quite well represented 
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by 



(21) 



where cq and ei are obtained by a fit to the simulation results (see Figj2|a)). 
Table |III| summarizes energy values obtained at fixed layer separation h = 1.05 for various 
dipole moments in the density range p = 0.3 — 0.7. For all densities considered the intralayer 
energy decreases with p and saturates near ^ ~ 2.5. The variation with density diminishes 
when the dipole moment is increased. The interlayer energy is much smaller than the intralayer 
contribution presenting, at all densities, a shallow minimum in the range ~ 1.75 — 2.0 where 
appreciable chaining of the particles sets in. 

C Pressure and surface stress 

Similar to the interlayer energy, the normal pressure at constant /i and p is quite well 
represented, as a function of /i, by 

However, as for a thermodynamical variable X generally 

/dX\ d<X > 
\l)h) ^ dh ' 

the fitting parameters /o and /i for the pressure do not relate directly to those for the energy. 



Nevertheless, the functional form of Eq.(22) obtained as the derivative of Eq.(21) provides quite 



good agreement between simulation results and Eq.(22) (see Fig|2|b)). 

As seen in Table |ITJ the surface stress, for p = 0.7, is fairly independent of h for p = 1 and 2. 
For p = 1, all the thermodynamic quantities, nlf"^^ and that contribute to fj through 



Eq.(B.15) and ([s]), are nearly constant. For p = 2, fj appears also to be insensitive to /i, but a 



small counterbalance between n^'^^ and 11^'^^ is observed as h increases from 1.01 to 1.15. As 
apparent from the one body orientational distribution functions, for h between 1.01 and 1.15 
and p = 2, the dipoles are on average less parallel to the layers than would be the case for 
larger h values. Thus, the attraction between particles in the same layer is slightly decreased in 
comparison to a monolayer ; this increases 11^'^^ and reduces Il'>^^\ since less contact between 
particles are observed in gi^trai^)- should note, though, that this effect is quite small (see 
Table |II|. 
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The values of fj, for fj, = 1, p = 0.7 and h > 2.00, given in table|ITj agree with the results obtained 
for the 2D pressure of the monolayer (see Tables I and II in ref . [Uj - as outlined in subsection 
3.2, the value of 11^ obtained from fj is twice the value of the pressure found in ref.[14!)- 
As shown previously, the 2D pressure of a monolayer of DHS may be related to the internal 
energy of the monolayer (see Eq.(21) in ref.[l3]). For the bilayer, we obtain almost exactly 
the same result, except for a factor 2 discussed before in subsection III.B. In Figjsj^a), we have 
represented —11^'^^ as a function of —f/™*'''^/^ j it appears that the dipolar contribution to the 
lateral pressure of the bilayer is very well represented by 

OTjintra jjintra 

li^^''^ = 3pkT ^^^^ = 3^^. (23) 

Thus, for p < 0.7 and p < 2.5, the equation of state is given by an equation similar to Eq.(21) 
of ref. [14 as 



2pkT pkT 2 iV 2pkT 

The variation of rj with dipole moment is shown in Fig. [sj^b) for h = 1.05 and various densities. 
fj can be approximated empirically by relations as 

flip, p) = -2pkT - 2U^^^\p, 0) + 5(ai; p, p) (25) 

where g{ai;p,p) is a function of the fitting parameter oi and IIj, (p, 0) obtained from the 
equation of state of hard disks (see, for instance, ref. [25]). Several functional forms for g, as for 
instance, gi{ai;p,p) = aip^/i^/(l + p'^), with oi ~ 2.7, or g2{a2;p,p) = a2/0^/U^/^, with 02 ~ 1.6 
were found to reproduce quite accurately the numerical results given in Table IV. 



D Structural properties 

Structural properties of the bilayer can be conveniently characterized by the coefficients 
/iiio, and /i220 of the expansion of the intra- and interlayer pair correlation functions 
^mtra(l)2) and hinteri^,'^) on a set of rotational invariants as described in subsection 3.3. Se- 
lected results for both intra- and interlayer correlation functions for h = 1.05 at densities p = 0.3 
and p = 0.7 are shown in Figsj4]-|6j The intralayer correlation functions for p = 1, reported in 
Fig. |4j agree very well with the correlation functions of the monolayer for the same p and for 
p = l (see Fig. 4 of ref.fH]). 
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The intralayer correlation functions present a succession of well defined peaks reflecting the for- 
mation of chains as also apparent from snapshots of configurations (FigsjTj^a) and[7|^b)). The 
peaks sharpen with increasing dipole moment indicating stronger bonding of the particles in the 
chains. The intralayer correlations appear to be quite insensitive to the layer separation and 
coincide within statistical error in the range h = 1.05 — 2.0. 

The interlayer correlation function gives information on the organization of particles in one layer 
relative to those in the other layer. Although the energy coupling between the layers is quite 
small one observes a strong correlation of the positional and orientational order of the particles 
in the two layers (at least for h < 2). Inspection of the interlayer distribution function gf^^^^ 
reveals, for dipole moments // > 2, a high probability of the particles to be on top of each other 
with opposite directions of the dipole moments {hj^^g^^ negative at s = 0). In addition, at dipole 
moments fj, > 2.25, peaks appear in g^^f^^ at s = (0.5 + n)a, {n = 0, 1,2...) at which hj^^^^ is 
positive giving evidence for configurations in which two chains in different layers are nearly on 
top of each other (possibly some lateral displacement) such that the chain axes of the two chains 
are displaced by half a HS diameter. In this case dipole moments point in the same direction. 
The effect is most pronounced at the lower density p = 0.3. 

The knowledge of hj^^^^{s) and hl^^^^{s) enables to recover intralayer and interlayer energies 
according to 

/ iQfrintra 1 

■112 ^(,) ds 



( fill 



N 

Pointer 

N 



O .10 



2 intra\ 



(26) 



2tt 



(s2 + /l2)3/2 



Similarly, the pressure tensor components are given by 



pidd) 



n 



(dd) 



-27r/iV2 



(s2 + /l2)5/2 



(27) 



oo 

S 



1"° 

Jo 



h}nter{s) ds 



10 - Jo {s^ + h^fl^ 

The quantities JJi-^^^"-^ jjmter^ Pzlf^ and H^'^^ computed with functions h\^^^[s) and /i, 



112 



''inter ('^)' 

can serve as a consistency check with the direct simulation results for energy and pressure using 



Ewald summations (Tables III and IV). Such a comparison is, however, conclusive only if the 
correlation functions decay to zero on the scale of the simulation box which was only fulfilled 
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at the lower fj, values (cf. Figs J4] - |6] for the correlation functions). For example at h = 1.05, 
p = 0.7 and /i = 1.0 one has ^C7^"*™/iV = -0.55, /3[7^"*^7iV = -0.16, Pjf ^ = -0.43 and 



_ 2Q in good agreement with the results of Tables III and IV, For h = 1.05, p = 0.7 and 



p = 2.0, integrating up to half the box length, one has /3C/™*™/A^ = -5.9, /5C/*"*''7A^ = -0.42, 
Pzz'^'^ = —1.40 and H^'^^ = —12.6 which compares favorably with the values of Tables III and 



Eqs. (26)-(27), show that we have the relation 11^'^^ = 3 f/*"-*'"'^/^ for /i — > oo ; this asymptotic 



behavior is in accordance with Eq.(23). However, it is surprising that Eq.(23) is verified with 



such accuracy even for /i = 1.05 (see subsection IV. C and Figjsj^a)). 



The values of /i^^*^ for s > 7 agree well with Eq.(18). For example, at p = 2.5 on has P2 ~ 0.42 



for both densities 0.3 and 0.7. This low value of P2 merely indicates some prevelant local nematic 
ordering but no global long range nematic ordering of the dipole moments. 
The characterization of the structural organization of the particles in the bilayer at high densities 
is subject to greater uncertainty due to system size dependence and convergence problems. To 
illustrate the difficulties we refer to snapshots of configurations at p = 0.9, p = 2 and h = 1.05 
taken at different "time" intervals during the MC evolution of the system shown in Figsjsj^a- 
d). The system, with 2 x 1600 particles, was started from two square lattices with random 
orientations of the dipole moments. Already after 500 cycles of trial moves small vortices have 
built up predominantly around particles with dipole moments oriented perpendicularly to the 
layers (Figjsf^a)). As sampling proceeds the vortices grow bigger and large patches develop 
within which particles arrange with local hexagonal order and parallel alignement of the dipole 
moments (Figj8][|b,c)), clearly an energetically favorable ordering. It remains somewhat unclear 
whether, for small system sizes, the p.b.c. can stabilize such a ferroelectric arrangement. Such 
a possibility was indeed observed for a smaller system size (2 X 576 particles) (cf. Figjsj^d)), 
and in one instance (/i = 1.005, p = 2) also for the 2 x 1600 system though an independent 
run of similar length (1 x 10^ cycles) at the same state point retained a vortex arrangement. In 
some cases, for the smaller 2 x 576 system, we also observed formation of stripes with opposite 
directions of the dipole moments. 

The structural behavior just described seems typical for dipole strength /i ~ 2 and not to depend 
much on layer separation in the range h = 1—2. For larger dipole moments the vortex structure 
appears to be more stable but, evidently, relaxation of the dipole moments is also slower. For 
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sure is that there are strong structural correlations between the layers. As for the lower densities, 
particles arrange preferentially to sit on top of each other with opposite directions of the dipole 
moments. 

Finally, in Figj9]we show the organization of dipole moments in a bilayer with h = 1.05 for close 
packed square and hexagonal lattices of the HS (disks). In both cases the HS in the two layers 
were taken to be on top of each other. On the square lattices {p = 1.0) the dipole moments 
in each layer align in parallel lines along the box edges with opposite directions of the dipole 
moments in neighboring lines (Figjoj^a)). A small tendency of microvortex formation is observed. 
These arrangements are typical of (monolayer) ground state configurations. For a square lattice 
of in-plane dipoles the ground state is continuously degenerated but thermal contributions can 
select configurations where rows or colums of parallel spins alternate [36]. In contrast, for the 2D 
triangular lattice with in-plane dipoles, the ground state of the infinite system is ferroelectric 
|37| |38] : in finite systems the dipolar ordering in the ground state may, however, depend on 
system size and aspect ratio of the lattice [39 . In the present finite temperature calculations 
(p = 1.15) we observe a ferroelectric phase with slight zigzag ordering of the dipole moments 
(Figjoj^b)). The influence on ordering of dipole strength, system size and use of p.b.c. has still to 
be investigated. It should be noted also that in our calculations the dipoles are not completely 
in-plane. As expected, for both lattices, dipole moments in different layers run in opposite 
directions. 

V SUMMARY AND CONCLUSION 

We have investigated by MC simulation the structural and thermodynamic properties of 
fully orientable dipolar hard spheres mobile in two parallel planar surfaces with particular em- 
phasis on the forces between the two layers. Interlayer correlations turn out to be quite small 
vanishing practically at layer separations of two HS diameters. The interlayer energy is attrac- 
tive for all states considered and the normal pressure is negative meaning that an external force 
must be supplied to keep the layers apart. Indeed isobaric MC simulations, allowing h to fluc- 
tuate, did not enable to find an equilibrium state; the system either collapsed (at low applied 
negative pressure) or the two layers drifted away (at larger pressures). The normal pressure is 
well described by a dependence at larger separations in agreement with a second order 
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perturbation theory of the interaction free energy of the surfaces in an infinite dielectric medium 
by Attard and Mitchell |33| I34j . Despite the weak interlayer energy there are strong correlations 
for the structural behavior of the particles in the two layers. Particles preferentially sit on top of 
each other with opposite orientations of the dipole moments. At densities of the order p ~ 0.9 
convergence of the MC sampling is slow and, moreover, finite size effects may affect the results. 
Although we believe that for large systems vortex formation is the preferred structure, arrange- 
ments with ferroelectric ordering or stripes with up and down orientations of the dipole moments 
were stabilized in the smaller systems, likely by the use of periodic boundary conditions. These 
problems clearly need a more detailed investigation. 

As an extension of the present work it would be of interest to consider the case where the media 
on either side of the layers have different dielectric constants, as would be the case, for instance, 
in a lipid bilayer model where the hydrocarbon tails and aquous regions are approximated by 
ideal dielectrics. Although the surface polarization arising from the dielectric discontinuities 
can in principle be taken into account through dielectric images |l0] few simulation results have 
been presented so far [41] . Such simulations could valuably add to the comprehension of the 
origin of the repulsive "hydration" forces measured in phospholipid bilayers at short distances 
|42] . Existing theoretical approaches based on continuum electrostatics |34| H3] seem to fail to 
predict correctly these repulsive forces. 
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APPENDIX A: EWALD SUMS FOR THE DIPOLAR 
ENERGY OF THE BILAYER 



The total dipolar energy of the bilayer computed with the Ewald method is written as 



U,, = Er + i^g^^ + E^g^^ + + i^G^o- 



Here Er is the short range (direct space) contribution to the energy given by 



(A.l) 



(A.2) 



with 



^, , erfc(Qr) 2a exp(— a^r^) 
B(r) = + 



IT 



(A.3) 



C(r) 



„erfc(ar) 2a/ n 3\exp(— a^r^) 

3 ^ + -^{2a^ + 



In Eq.(A.2) it is assumed that the parameter a is sufficiently large to restrict interactions to 



the basic simulation cell. The energy Er can, in turn, be separated into an intralayer E^'^'^"' 
and an interlayer E'^*^'" contribution. The four last terms in Eq.(A.l) are the reciprocal space 



contributions. Each of the terms is again separated into intralayer and interlayer contributions. 
They are split into three contributions: Eq^^ involves only coupling between the normal com- 



(2) 

ponents of dipole moments, Eq_^^ coupling between in-plane and normal components of dipoles 
(3) 

and Eq_^^ in-plane coupling. Contributions to the interlayer energy are given by 

7 I(a,G;/i) X SRe[( ^ /if exp(iG- Si))( ^ exp(-iG • s^)) 



,(1, inter) 



ieLi 



E 



(2,mter) _ TT 



G^O 



A 



G^o 



ieLi 



IJ-i exp{iG ■ Si))[ Y, {f^j • G) exp(- iG ■ sj))] 



E. 



{3,inter) 

G^o 



J Y K(a, G; h) x 3f?e[( Y ' G) exp(iG • s,)){ Y if^j ' G) exp(-iG • s,))] 



G^o 



(A.4) 
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where ?R.e[z] and 9m[z] are the real and imaginary parts of the complex number z, respectively. 
G = 2tt(j^, jj-)-, {uxjUy integers) is a two-dimensional vector in recriprocal lattice and G =|| G \\. 
The functions I(a, G; h), J(a, G; h) and K(a, G; h) are given by 

I(a,G; h) 



A 

" exp ( - - a^h^) - G'^K{a, G; h) 



G G 

3(a,G;h) = exp(G/i)erfc( \- ah) — exp(— G/i)erfc( ah] 

2a 2a 



(A.5) 



1 r G G y 

K{a,G;h) = — exp(G/i)erfc( h a/i) + exp(— G/i)erfc( ah 

G 2a 2a 

The constant term is 



Contributions to intralayer the energy are given by 



(A.6) 



E 



(l,intra) 

G^o 



D(«, G) [ I ^ //f exp(iG • Si) |2 + I exp(iG • Sj) \^ ] 



„{2,intra) „ 



with 



(3, intra) 

G^o 



H(a, G) [ I ^ (m, • G) exp(iG • s,) I' + I E (z^^' " ^) ^^P(^^ " 

G^o 



D(a,G) 
H(a,G) 



2a 



TT 



exp(-G74a^) - G erfc(G/2a) 



(A.7) 



erfc(G/2a) 
G 



and the constant is 



^(intra) ^ 2a^ ^ ^^^^2 ^ ^ ^ ^^^^2] 2a ^ __2 



'G=o 



^7. 



(A.9) 



Due to the 2d character of G it is easily seen from the corresponding term in Eq.(A.4) (interlayer 

p(2, intra) 



contribution) that Eq_^^ must vanish. 
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APPENDIX B: THE MICROSCOPIC STRESS TENSOR OF 

THE BILAYER 



In this Appendix, we derive the microscopic stress tensor for the bilayer system from its equations 
of motion, in a way similar to the one of ref.[26](a) for inhomogeneous fluids. The microscopic 
stress tensor of the bilayer is split into normal and lateral or components as 



(J = (7t + <7 Af 



XX ^xy 



\ 



V 



+ 











^xz 
O'yz 



\ 



(B.l) 



y (^xz ^yz ^zz J 





The Lagrangian function of the bilayer system, with the constraints Zi = Hi, for i G Li, and 
Zi = H2, for i G L2 IS given by 



miH2 



ieLi 



(B.2) 



where $ is the pair potential energy due to interactions between particles and ^ext represents 
the action of any external fields. In the above equation. Hi and H2 are collective variables 
associated with the z-coordinate of the layers. From the lagrangian of the system, we obtain the 
equations of motion for the particles in the layer Li and the collective variable Hi : 



- Vi$(Sij,0) - Y V,$(Sij,i72 - Hi) - Vi<^ext{Si,Hi] 



and similar equations for the layer L2. m denotes the mass of the particles. 



(B.3) 
(B.4) 



The momentum density for the bilayer system can be written as 
J{s,z,t) = JT{s,z,t) + JN{s,z,t)e;, 



m 



6{z - Hi) ^ Si 6{s - Si) +m 6{z - H2) ^ Sj 6{s - Si 



(B.5) 



+mHi 5{z - Hi) ^ 5{s - Si)ez + mH2 S{z - H2) ^ S{s - Si)ez 

ieLi i&L2 
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where 6{x) is the Dirac distribution. From the time derivative of the momentum density, we 
obtain easily [26] the kinetic contribution to the lateral component of the stress tensor as 

a^p{s, z, t) = -m 6{z - Hi) s^''si^ 6{s - Si) - m 6iz - H2) ^ Si°s/ 6{s - s,) (B.6) 
with a,P = X, y. The kinetic contribution to the normal component is obtained similarly as 



cjf,(s, z, t) = -mHi S{z - Hi) ^i" ^(^ " ^i) " S{z - H2) ^ Si" 5{s - Si 



a^^{s,z,t) = —mHi' 5{z — Hi) 6{s — Si) — mH2^ S{z — H2) 6{s — Si 



The configurational contributions to the stress tensor, follow from Eq.(B.3) 



(B.7) 



S{z - Hi 



+ 



(B.8) 



+ ^ E E - H2) [ df 6{s - I) 



5{z - H2) 



with a = x,y and Cij a contour joining Sj to Sj in the plane perpendicular to the z direction. 



Eqs.(B.8) and (B.6) allow to fully determine the lateral component of the stress tensor of the 



bilayer. The integrals in Eq.(B.8) can be evaluated by using the parametrization proposed by 



Irving and Kirkwood [26j(b), namely 

^ ^ Vfc^(s,„0) / df5{s-l) = Y, E 4^?'^{^ij,0) [' dX6{s-Xsj-{l-X)s,) 



(B.9) 
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and 



^ y?Hs,j,H2 - Hi) I' df 5{s - I) 



E 



^2 - ^i) /' d\ 6{s - Xs, - (1 - \)s,) (B.IO) 



y. y. 4^?HsiJ'Hi - H2) C d\ 5{s - Xsj - (1 - X)si 



Eqs.(B.6) and (B.8) show that can be written in the form (a, P = x,y) 



a^p{s,z,t) = ri^'^{s,t) 6{z-Hi) + T^^^{s,t) 6{z - H2) (B.ll) 

One should note that, if z 7^ Hi and z 7^ H2 then (Tq,^(s, z, t) = 0. 

In accord with sohd surface physics we define the surface stress tensor as 

Vap{s,t) = j ao,p{s,z,t)dz = T^^^{s,t)+T^^^{s,t) (B.12) 

If one adopts the two-component monolayer picture discussed in the main text, then each contri- 
bution rj^i^ and rj^^ correspond respectively to partial contribution of each species to the surface 
stress tensor. 

From the surface stress tensor we define the lateral component of the pressure tensor of the 
bilayer as the ensemble average of the surface stress tensor as 



It follows that 



^ / dsijo.p{s,t)). (B.13) 
^ JL1UL2 ' 



7E E4V"^(^^^'^2-^l) 



A 

The average lateral pressure 11^ and the surface stress 77 are then given by 



UT = -iIl^x+Ilyy) = -f] (B.15) 
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The configurational contribution to the normal component Ozz ahows to obtain the force acting 
on the layers. From the equations of motion of H\ and 7^2; we obtain 

(B.16) 

Thus, the total force F2_,^ acting on layer Li due to the particles in layer L2 is given by 



and, obviously, we have 



>2 



ds— z = H2,t) = -Fi_^ (B.18) 



The average force by unit area is 



The equation (B.19) for P/v is in full agreement with the derivation of the normal pressure 



derived for similar systems in refs.j26l [271 1 
If the z-coordinates of the layers are fixed, as is the case in most of the computations in the 
present work, an external field compensates exactly the microscopic forces. In this case we have 
Hi = —H2 = h/2, Hi = H2 = and Hi = H2 = and the external forces are given by 

^-M = E ^^-*(^- ^) = -^2-1 (B-20) 

and 
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APPENDIX C: RECIPROCAL SPACE CONTRIBUTIONS 
TO THE PRESSURE TENSOR AND FORCES 



The general formulae for the components of the stress tensor in terms of the interaction 
potential are given in section 2. In this appendix, we give explicit expressions for the reciprocal 
space contribution in an Ewald sum of the stress tensor components. They can be obtained 
directly from the results of Appendix A or from the general derivation given by Heyes |17] for 
quasi- two dimensional systems. 

The short ranged contributions are easily obtained from Eqs.(A.2-3). 

From Eq.Q and with notations of Appendix A, we have, for the bilayer system, 

nf '^^ = ^ • VsAE^G^^^ + E^G^^)) (C.l) 

i 

The intralayer contributions to the lateral components of the stress tensor are given by 



XKsm 



[{^(.G- Si)^il exp(zG • Si)){ ^if exp(iG • Si)) (c.3) 



+ ( J^(G-Si)/ifexp(iG-Si))( Y fil exp{iG ■ Si)) 



^..•v.,^g;7'^) = o (C.4) 



2tt 



[^iG- s,)ifi, ■ G) expiiG ■ s,)){ ^ (/x, • G) exp(iG • s,)) 



; 5^ (G • s,)i^Li . G) exp(iG ■s,)){ J] (/x, • G) exp(iG ■ s,))] 

(C.5) 
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with functions D and H as defined in Eq.(A.8) 
Interlayer contributions are given by 



4( ^^ieM^G • Si)) ( lM'j{G ■ sj)exp(-iG • s^)) (c.6) 



^ /if (G • Si)exp(zG • Si)) ( ^ /x^exp(-iG • Sj 



x3f?e[( /xf(G • Si)exp(iG • Si))( ■ G)exp(-iG • s,)) 

ieii i£-J^2 

-( X] /^fexp(iG • Si))( X (/Xj. • G)(G • Sj)exp(-iG • s^)) (c.7) 

«6Li ieL2 

+ ( ^(/Xi • G)(G • Si)exp(iG • Si))( ^ //|exp(-iG • s^)) 

«6Li ie-L2 



-( X(/x,, • G)exp(iG • Si))( J] /.|(G • s,)exp(-zG • s,)) 
vr 



.[( X(/x, • G)(G • Si)exp(iG • s,))( (/^, • G)exp(-iG • s,)) 



-( X (/X, • G)exp(iG • s,)) ( (/X, • G)(G • s,)exp(-iG • s,)) 



(C.8) 



with functions /, J and ii' defined in Eq.(A.5). 
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The contributions to the normal component of the stress tensor are given by 



^ inter) 



z=h 



A 



^ [G2j(a,G;/i) + ^Q(a,G;/i)] 



:5Re[( ^ /if exp{iG ■ Sj)) ( ^ exp(— iG • s 



d „{2,inter) 



z=h 



TT 

A 



2a 



^ [G^K{a, G; h) - ^P(a, G; /i) 



G^o 



exp(iG • Si)) ( {fij ■ G) exp(-iG • Sj))] 



^ 771(3, inter) 

d~z G^o 



z=h 



TT 

A 



G^o 



x^e[{ J2ifJ'i-G)exp{iG-Si)){ ^ (/i^ • G) exp(-iG ■ 
The function Q{a, G; h) is obtained from the derivative of J, i. e. 

q2 

Q(q, G] h) = 2exp ( k) exp(— Q^/i^). 

4a^ 



Finally, 



dz G=o 



z=h 



G=o 



with e'q^^'^ given by Eq.(A.6) 
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List of Tables 



Table |I] : Definitions of tlie projections of intralayer and interlayer correlation func- 
tions computed in the present work. 



Table II : Average energies and pressures for the bilayer system for p = 0.7 and 
several values of h. The numbers in brackets give the accuracy on the last digit of the 
averages, a is the width of the one-body orientational distribution obtained by fitting the 
MC histograms. (3Udd/N, /5?7*"*'""/A^ and PU^"*^^ /N denote, respectively, the averages of 
total, intralayer and interlayer dipolar energies. Pif'^'' is the average normal force by unit 
area as defined by Eq.(j7|. H^'^^ is the average of the dipolar contribution to the lateral 
pressure computed with Eq.(|4) and 11^'^^ is the hard sphere contribution computed from 
the contact value of the pair distribution function Eq.(6). fj = —2pkT — 2Il'>^^^ — 11^'^^ is 



the surface stress as defined in Eq.(B.15) 



Table III : Average energies for the bilayer system for several values of p and p for 
h = 1.05. Notations are the same as in Table HTl 



Table IV : Average pressures for the bilayer system for several values of p and p for 
h = 1.05. Notations are the same as in Table HTl 
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List of Figures 



Figure [TJ Orientational distribution functions of dipolar moments in a bilayer of 



dipolar hard spheres. Symbols denote MC data and solid lines are fits using Eq.(19). (a) 
Results at selected values of p and fi aX h = 1.05. (b) Variation with layer separation h 
for p = 0.7 and p = 2.00. 

Figure [2| Average energies (a) and normal pressures (b) as a function of h for p = 0.7 
and /i = 1 and 2. The symbols denote MC data and the lines are fits to the data us- 



ing Eqs.(21) and (22), respectively. The fitting parameters for p=l are eo = 0.16 ± 0.01, 
ei = 0.045±0.002 and /o = 0.46±0.01, /i = 0.13±0.01. For p=2 they are eo = 0.31±0.01, 
ei = 0, 28 ± 0.01 and /o = 0.68 ± 0.03, /i = 1.3 ± 0.1. 

Figure |3| (a) Lateral pressure as a function of the intralayer energy per unit area. 
Symbols are data from Tables III and IV for densities p = 0.3 — 0.7, dipole strengths 



p = 1.0 — 2.5 and h = 1.05 ; the straight line is given by Eq.(23). (b) Surface stress 



as function of dipole strength for p = 0.3 — 0.7 and h = 1.05. Symbols are data from 



Table IV and lines are given by Eq.(25) with gi{ai; p, p) = aip^ p'^ / {1 + p^) (ai ~ 2.7) and 



n^*^^ (p, 0) given by the equation of state of hard disks ref . 



Figure [4| Intralayer angle averaged pair distribution function (a) g^ntrai.^) angu- 
lar projections (b) h^l^.^{s) for a bilayer of dipolar hard spheres at p = 0.7, h = 1.05 for 
p = 1.0 (black) and p = 2.0 (red -grey hue). 

Figure [5| Interlayer angle averaged pair distribution function ginteri^) angular 
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projections h^^^j.{s) of the pair distribution functions gmteri^'^) for the DHS bilayer at 
p = 0.3 and h = 1.05 for several values of p. (a) ^"Xrl^) ; (b) hj^teAs) ; (c) hj^l^s), (d) 

"' inter \^ )■ 

Figure |6j Same as Fig. [5] but for p = 0.7. 

Figure [7[ Snapshots of bilayer configurations of particles at /i = 2.0 ((a),(b)) and 
p = 2.50 ((c), (d)) for h = 1.05 ; snapshots (a) and (c) are for p = 0.3 (A^ = 1058); 
snapshots (b) and (d) for p = 0.7 (A^ = 1024). Particles in different layers are represented 
by different colours. The HS cores are represented by circles of diameter a = 1 and the 
directions of dipole moments by arrows. 

Figure |8} Bilayer configurations of the 2 x 1600 particle system at p = 0.9, p = 2 
and h = 1.05 at different intervals of the MC simulation; (a) snapshot after 500 cycles, 
(b) 0.26 X 10^ cycles, (c) 1.75 x 10^ cycles, (d) resuh for 2 x 576 particles after 2.6 x 10*^ 
cycles. For clarity only the particle arrangements in one layer are shown in (a)-(c). The 
arrows denote the projections of the dipole moments on the layer plane. Thus dipoles 
perpendicular to the layer appear as dots. 

Figure [9| Snapshots of bilayer configurations of particles at close packing, (a) square 
lattice {p=l,p = 2,h = 1.05, N=3200); (b) hexagonal lattice (p = 1.15, p = 2, h = 1.05, 
N=2400). The particles in the two layers are on top of each other. The arrows denote the 
projections of the dipole moments on the layer plane. The two layers are shown separately. 
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TABLE I: Definitions of the projections of intralayer and interlayer correlation functions 
computed in the present work. 
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TABLE II: Average energies and pressures for the bilayer system for p = 0.7 and several 
values of h. The numbers in brackets give the accuracy on the last digit of the averages. 
a is the width of the one-body orientational distribution obtained by fitting the MC 
histograms. pUdd/N, /5?7*"*™/A^ and PW'^^'^^/N denote, respectively, the averages of total, 
intralayer and interlayer dipolar energies. Pif) is the average normal force by unit 
defined by Eq.([7|. II^'^^ is the average of the dipolar contribution to the lateral pressure 

AHS) 



computed with Eq.(|4) and II^^ is the hard sphere contribution computed from the 
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contact value of the pair distribution function Eq.(6|. rj = —2pkT — 2Il'>^^^ — II^'''' is the 



surface stress as defined in Eq.(B.15). 
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TABLE III: Average energies for the bilayer system for several values of p and /i for 
h = 1.05. Notations are the same as in Table [TTl 
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TABLE IV: Average pressures for the bilayer system for several values of p and [i 
h = 1.05. Notations are the same as in Table [TTl 
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FIG. 1: Orientational distribution functions of dipolar moments in a bilayer of dipolar 



hard spheres. Symbols denote MC data and solid lines are fits using Eq.(19). (a) Results 
at selected values of p and fi at h = 1.05. (b) Variation with layer separation /i for p = 0.7 
and fi = 2.00. 
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FIG. 2: Average energies (a) and normal pressures (b) as a function of /i for p = 0.7 
and /i = 1 and 2. The symbols denote MC data and the lines are fits to the data using 
Eqs.(21) and (22), respectively. The fitting parameters for /i=l are cq = 0.16 ± 0.01, 
ei = 0.045±0.002 and /o = 0.46±0.01, /i = 0.13±0.01. For fi=2 they are Cq = 0.31±0.01, 
ei = 0, 28 ± 0.01 and /o = 0.68 ± 0.03, fi = 1.3 ± 0.1. 
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FIG. 3: (a) Lateral pressure as a function of the intralayer energy per unit area. Symbols 
are data from Tables III and IV for densities p = 0.3 — 0.7, dipole strengths fi = 1.0 — 2.5 



and h = 1.05 ; the straight line is given by Eq.(23). (b) Surface stress as function of 
dipole strength for p = 0.3 — 0.7 and h = 1.05. Symbols are data from Table IV and lines 
are given by Eq.(25) with gi{ai;p,p) = aifB^'^/{l + p"^) (ai ~ 2.7) and ll^^\p,0) given 
by the equation of state of hard disks ref. 
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FIG. 4: Intralayer angle averaged pair distribution function (a) Qintrai^) angular 
projections (b) /iSa(s) for a bilayer of dipolar hard spheres at p = 0.7, h = 1.05 for 
/i = 1.0 (black) and /i = 2.0 (red -grey hue). 
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FIG. 5: Interlayer angle averaged pair distribution function ginteri^) angular projec- 
tions /if^terl-^) of the pair distribution functions (yfj„ier(12) for the DHS bilayer at p = 0.3 
and h = 1.05 for several values of /i. (a) g^Z,{s) ; (b) hj^J^^.{s) ; (c) hj^l^{s), (d) h^Z^.{s). 
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(a) (b) 




FIG. 7: Snapshots of bilayer configurations of particles at /i = 2.0 ((a),(b)) and /i = 2.50 
((c), (d)) for h = 1.05 ; snapshots (a) and (c) are for p = 0.3 {N = 1058); snapshots (b) 
and (d) for p = 0.7 (A^ = 1024). Particles in different layers are represented by different 
colours. The HS cores are represented by circles of diameter a = 1 and the directions of 
dipole moments by arrows. 
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FIG. 8: Bilayer configurations of the 2 x 1600 particle system at p = 0.9, 11 — 2 and 
/i = 1.05 at different intervals of the MC simulation; (a) snapshot after 500 cycles, (b) 
0.26 X 10^ cycles, (c) 1.75 x 10^ cycles, (d) result for 2 x 576 particles after 2.6 x 10^ 
cycles. For clarity only the particle arrangements in one layer are shown in (a)-(c). The 
arrows denote the projections of the dipole moments on the layer plane. Thus dipoles 
perpendicular to the layer appear as dots. 
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FIG. 9: Snapshots of bilayer configurations of particles at close packing, (a) square lattice 
(p = 1, = 2, /i = 1.05, N=3200); (b) hexagonal lattice (p = 1.15, = 2, /i = 1.05, 
N=2400). The particles in the two layers are on top of each other. The arrows denote 
the projections of the dipole moments on the layer plane. (The two layers are shown 

separately) . 
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